function [out, Ct, Pt, Dt, Qt, Gt] = ll(xx, T_f, SET, data)
%
% Calculates the log-likelihood for Gaussian sequence given a mean and
% variance sequence.
%
% State space representation:
% 
% S_t = C_t + P_t S_{t-1} + D_t n_t : State-state equation
% X_t = H S_t + v_t : observation equation

n_   = SET.variable.n_ ; 
nobs = SET.nobs ;

%% Re-solve model with xx parameters

if or(sum(xx<0)>0,sum(xx(1:5)>1)>0)
    out=-1e12;
    return
end

SET.M_.params(1:SET.EST.nparam_est) = xx ;

SC.t_f = T_f ;

Tbstar = SET.tv_str_chg.Tbstar ;

if Tbstar>0
    SET.tv_str_chg.cur_t  = 1 ; 
    
    mats = tvmats(SET,SC) ;
    
    if mats.unq==0 
        out=-1e12;
        disp('not unq') ;
        return
    end
    
    Qt = mats.Qt ; 
    Gt = mats.Gt ; 

    Qt = Qt(:,:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date));
    Gt = Gt(:,:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date));
else
	dyn_in.M_  = SET.M_ ;
    dyn_in.oo_ = SET.oo_ ;
    dyn_in.options_ = SET.options_ ;
    dyn_in.solve = 1 ;
    dyn_in.linearize_around_diff_y = 0 ;

    mats = resolve_(SET,dyn_in) ;
    
    if isempty(mats.mats.Q)
        out=-1e12 ;
        disp('not unq') ;
        return
    end
    
    Qt = repmat(mats.mats.Q,[1 1 SET.ss]) ;
    Gt = repmat(mats.mats.G,[1 1 SET.ss]) ;
    
    zlb_t = SET.zlb_t ;

    SET.mat_i_f_zlb = mats.mat_i_f_zlb ;
    SET.mat_init = mats.mat_init ;
    SET.mats.Q = mats.mats.Q ;

    if max(T_f)>0 
        zlb_mats = rfmats(SET, 1, max(T_f)) ; % Solve for reduced form ZLB mats
        tmp = max(T_f)-T_f+1 ; % Index
        Qt(:,:,zlb_t>0) = zlb_mats.Qhat(:,:,tmp(zlb_t>0)); % Insert ZLB in Qt
        Gt(:,:,zlb_t>0) = zlb_mats.Ghat(:,:,tmp(zlb_t>0)); % Insert ZLB in Gt
    end
end    
    
%% Extract matrices for Kalman filter

Ct = Qt(1:n_-1,end,:) ;
Pt = Qt(1:n_-1,1:n_-1,:) ;
Dt = Gt(1:n_-1,:,:) ;

H = zeros(nobs,n_-1) ; 

for i=1:nobs
    H(i,SET.EST.obs_(i)) = 1 ;
end

[U, St, flag] = kf_(SET, Ct, Pt, Dt, H, data, SC) ;

if flag ; out=-1e12 ; return ; end

out = 0;

for t=1:SET.ss
    hidx = SET.EST.hidx ;
    if SC.t_f(t)>0 ; hidx = SET.EST.hidx_no_i ; end
    
    out = out - (nobs/2) * log(2*pi) - ...
        (1/2) * log(det(St(hidx,hidx,t))) - ...
        (1/2) * U(hidx,t)' / St(hidx,hidx,t) * U(hidx,t) ;
end

if abs(imag(out))>0
    out=-1e12;
    disp('imag ll') ;
end
